
%outputs the ellipse parameterization and the confidence interval
function [r,ci] = refineHoughEllipse(po,ridgeinds,imsz,bounds)

if length(imsz) == 3
    nf = imsz(3);
else
    nf = 1;
end
r = nan(nf,5);
ci = nan(nf,5,2);
%get lower and upper bounds for ellipse size based on the range of initial
%parameters:
if nargin < 4 || isempty(bounds)
    lb = zeros(1,5); ub = lb;
    ellipax1 = po(:,3,:);
    ellipax2 = po(:,4,:);
    ellipax1 = nanmean(ellipax1(:));
    ellipax2 = nanmean(ellipax2(:));
    if max(ellipax1,ellipax2)/min(ellipax1,ellipax2) > 2.5
        temp = po(:,3,:);
    else
        temp = po(:,3:4,:);
    end
    lb(1) = min(temp(:)*.8);
    if lb(1)<0
        lb(1) = 0;
    end
    ub(1) = max(temp(:)*1.25);

    lb(2) = .85; ub(2) = 1.17;       %more conservative...
    
    
    temp = po(:,1,:);
    lb(3) = min(temp(:))*.8; ub(3) = max(temp(:))*1.25;
    temp = min(po(:,2,:));
    lb(4) = min(temp(:))*.8; ub(4) = max(temp(:))*1.25;
    lb(5) = -pi; ub(5) = pi;
    
    bounds = [lb;ub];
else
    lb = zeros(1,5); ub = lb;
    lb(3) = min(po(:,1))*.7; ub(3) = max(po(:,1))*1.35;
    lb(4) = min(po(:,2))*.7; ub(4) = max(po(:,2))*1.35;
    bounds(:,3:4) = [lb(3:4);ub(3:4)];
    bounds(1,2) = .85;
    bounds(2,2) = 1.17;

end

for n = 1:nf
    ind = ridgeinds(:,n);ind = ind(~isnan(ind));
    [yridge,xridge] = ind2sub(imsz(1:2),ind);
    eo = po(n,[3,4,1,2,5],1);eo(5) = deg2rad(eo(5));
    
    % check if you're inputting minor and major axes or one axis + aspect
    % ratio (this is sort of a patch:
    if eo(2) > 5
        eo(2) = eo(2)/eo(1);
    end
    if eo(1)<bounds(1,1) || eo(1) > bounds(2,1)
        eo(1) = (bounds(2,1)-bounds(1,1))/2 + bounds(1,1);
    end
    if eo(2)<bounds(1,2) || eo(2) > bounds(2,2)
        eo(2) = (bounds(2,2)-bounds(1,2))/2 + bounds(1,2);
    end
    if eo(5)<bounds(1,5) || eo(5) > bounds(2,5)
        eo(5) = (bounds(2,5)-bounds(1,5))/2 + bounds(1,5);
    end
    
    
    try
        [params,curci] = optimizeEllipse(eo, [yridge,xridge],bounds);
        r(n,:) = params;
        ci(n,:,1) = curci(:,1); ci(n,:,2) = curci(:,2);
        
    catch
        r(n,:) = nan;
        disp('unable to refine')
    end
    disp(n/nf)
end

